home *** CD-ROM | disk | FTP | other *** search
/ Celestin Apprentice 5 / Apprentice-Release5.iso / Source Code / Libraries / LinAlg 3.1 / LinAlg / vsvd.cc < prev    next >
Encoding:
C/C++ Source or Header  |  1995-12-21  |  4.0 KB  |  141 lines  |  [TEXT/ALFA]
  1. // This may look like C code, but it is really -*- C++ -*-
  2. /*
  3.  ************************************************************************
  4.  *
  5.  *    Verify Singular Vector Decomposition of a Rectangular Matrix
  6.  *
  7.  * $Id$
  8.  *
  9.  ************************************************************************
  10.  */
  11.  
  12. #include "svd.h"
  13. #include <iostream.h>
  14.  
  15.                 // SVD-decompose matrix A and test we can
  16.                 // compose it back
  17. static void test_svd_expansion(const Matrix& A)
  18. {
  19.   cout << "\n\nSVD-decompose matrix A and check if we can compose it back\n"
  20.        << endl;
  21.   A.print("original matrix");
  22.   SVD svd(A);
  23.   svd.q_U().print("left factor U");
  24.   svd.q_sig().print("Vector of Singular values");
  25.   svd.q_V().print("right factor V");
  26.   
  27.   {
  28.     cout << "\tchecking that U is orthogonal indeed, i.e., U'U=E and UU'=E"
  29.          << endl;
  30.     Matrix E(Matrix::Unit,svd.q_U());
  31.     Matrix ut(Matrix::Transposed,svd.q_U());
  32.     Matrix utu(svd.q_U(),Matrix::TransposeMult,svd.q_U());
  33.     verify_matrix_identity(utu,E);
  34.     Matrix uut(svd.q_U(),Matrix::Mult,ut);
  35.     verify_matrix_identity(uut,E);
  36.   }
  37.   
  38.   {
  39.     cout << "\tchecking that V is orthogonal indeed, i.e., V'V=E and VV'=E"
  40.          << endl;
  41.     Matrix E(Matrix::Unit,svd.q_V());
  42.     Matrix vtv(svd.q_V(),Matrix::TransposeMult,svd.q_V());
  43.     verify_matrix_identity(vtv,E);
  44.     Matrix vt(Matrix::Transposed,svd.q_V());
  45.     Matrix vvt(svd.q_V(),Matrix::Mult,vt);
  46.     verify_matrix_identity(vvt,E);
  47.   }
  48.   
  49.   {
  50.     cout << "\tchecking that U*Sig*V' is indeed A" << endl;
  51.     Matrix S(Matrix::Zero,A);
  52.     (MatrixDiag)S = svd.q_sig();
  53.     Matrix vt(Matrix::Transposed,svd.q_V());
  54.     S *= vt;
  55.     Matrix Acomp(svd.q_U(),Matrix::Mult,S);
  56.     compare(A,Acomp,"Original A and composed USigV'");
  57.   }
  58.  
  59.   cout << "\nDone" << endl;
  60. }
  61.  
  62.                 // Make a matrix from an array
  63.                 // (read it row-by-row)
  64. class MakeMatrix : public LazyMatrix
  65. {
  66.   const REAL * array;
  67.   const int no_elems;
  68.   void fill_in(Matrix& m) const;
  69. public:
  70.   MakeMatrix(const int nrows, const int ncols,
  71.            const REAL * _array, const int _no_elems)
  72.     : LazyMatrix(nrows,ncols), array(_array), no_elems(_no_elems) {}
  73.   MakeMatrix(const int row_lwb, const int row_upb,
  74.            const int col_lwb, const int col_upb,
  75.            const REAL * _array, const int _no_elems)
  76.     : LazyMatrix(row_lwb,row_upb,col_lwb,col_upb),
  77.       array(_array), no_elems(_no_elems) {}
  78. };
  79.  
  80. void MakeMatrix::fill_in(Matrix& m) const
  81. {
  82.   assert( m.q_nrows() * m.q_ncols() == no_elems );
  83.   register const REAL * ap = array;
  84.   for(register int i=m.q_row_lwb(); i<=m.q_row_upb(); i++)
  85.     for(register int j=m.q_col_lwb(); j<=m.q_col_upb(); j++)
  86.        m(i,j) = *ap++;
  87. }
  88.  
  89. static void test1(void)
  90. {
  91.   cout << "\nRotated by PI/2 Matrix Diag(1,4,9)\n" << endl;
  92.   
  93.   REAL array[] = {0,-4,0,  1,0,0,  0,0,9 };
  94.   test_svd_expansion(MakeMatrix(3,3,array,sizeof(array)/sizeof(array[0])));
  95. }
  96.  
  97. static void test2(void)
  98. {
  99.   cout << "\nExample from the Forsythe, Malcolm, Moler's book\n" << endl;
  100.   
  101.   REAL array[] = 
  102.        { 1,6,11, 2,7,12, 3,8,13, 4,9,14, 5,10,15};
  103.   test_svd_expansion(MakeMatrix(5,3,array,sizeof(array)/sizeof(array[0])));
  104. }
  105.  
  106. static void test3(void)
  107. {
  108.   cout << "\nExample from the Wilkinson, Reinsch's book\n" <<
  109.           "Singular numbers are 0, 19.5959, 20, 0, 35.3270\n" << endl;
  110.   
  111.   REAL array[] = 
  112.       { 22, 10,  2,   3,  7,    14,  7, 10,  0,  8,
  113.         -1, 13, -1, -11,  3,    -3, -2, 13, -2,  4,
  114.          9,  8,  1,  -2,  4,     9,  1, -7,  5, -1,
  115.          2, -6,  6,   5,  1,     4,  5,  0, -2,  2 };
  116.   test_svd_expansion(MakeMatrix(8,5,array,sizeof(array)/sizeof(array[0])));
  117. }
  118.  
  119. static void test4(void)
  120. {
  121.   cout << "\nExample from the Wilkinson, Reinsch's book\n" <<
  122.           "Ordered singular numbers are Sig[21-k] = sqrt(k*(k-1))\n" << endl;
  123.   Matrix A(21,20);
  124.   struct fill_in : public ElementAction
  125.   {
  126.     void operation(REAL& element) { element = j>i ? 0 : i==j ? 21-j : -1; }
  127.   };
  128.   A.apply(fill_in());
  129.   test_svd_expansion(A);
  130. }
  131.  
  132. main(void)
  133. {
  134.  cout << "\n\nTesting Singular Value Decompositions of rectangular matrices"
  135.       << endl;
  136.  test1();
  137.  test2();
  138.  test3();
  139.  test4();
  140. }
  141.